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BOUNDARY-INTEGRAL METHODS IN ELASTICITY AND PLASTICITY 

by Alexander Mendelson 
Lewis Research Center 

SUMMARY 

Recently developed methods that use boundary -integral equations applied to elastic 
and elastoplastic boundary value problems are reviewed. Direct, indirect, and semi- 
direct methods using potential functions, stress functions, and displacement functions 
are described. Examples of the use of these methods for torsion problems, plane prob- 
lems, and three-dimensional problems are given. It is concluded that the boundary - 
integral methods represent a powerful tool for the solution of elastic and elastoplastic 
problems. 


INTRODUCTION 

Methods of analysis in solid mechanics, as in most other scientific and engineering 
fields, have been revolutionized by the advent of the modern digital computer. Thus, 
purely numerical methods, such as finite differences, and analytical methods, whose 
final stages require numerical procedures such as complex variable methods, awaited 
the modern computer for their practical implementation to all but the simplest problems. 

Although many general methods have been and are being used for solving problems 
in elasticity and plasticity, the currently most popular ones probably are (1) finite differ- 
ences, (2) finite elements, and (3) complex variables. To this list we may now add the 
boundary -integral methods. It is the boundary -integral methods that are discussed in 
this report. 

First, the question might be asked, why is it necessary to get involved in yet another 
method for solving the same types of problems ? What advantages do these methods have 
over the commonly used methods previously listed? The apparent advantages are as 
follows. 

(1) The need for conformal mapping is obviated. 

(2) Uniform or mixed boundary value problems are handled with equal ease. 



(3) Stresses and displacements are obtained directly without need for numerical 
differentiation. 

(4) No special treatment is needed for multiply connected regions. 

(5) The internal stresses and displacements are obtained only where and when 
needed. 

(6) The method can be extended directly to three-dimensional problems. 

(7) Nodal points are needed only on the boundary instead of throughout the interior 
(as required by finite -difference or finite -element methods). 

The last point, which is probably the most important one, is illustrated in figure 1. 
For the finite -difference or finite -element methods, the whole region must be covered by 
a grid, producing a large number of nodal points and corresponding unknowns. Thus a 
large number of simultaneous equations must be solved. In the present methods, as will 
shortly be seen, nodal points are taken only on the boundary, resulting in a much smaller 
number of unknowns, the dimensions of the problem being effectively reduced by one. 

On the other hand, the resulting matrices are full, whereas in the finite -element 
methods, for example, the matrices are usually sparse and can be more or less banded. 

The boundary integral methods themselves can be classified into three groups: 
indirect, semidirect, and direct. 

The indirect methods formulate the problem in terms of unknown density functions, 
which have no physical significance. But once these density functions are determined, 
the displacements and/or stresses can be directly computed. 

The semidirect methods formulate the problem in terms of unknown functions that 
are more familiar, such as the various stress functions. The stresses are then deter- 
mined by simple differentiation. 

The direct methods formulate the problem in terms of the direct physical quantities 
such as the displacements. 

All three methods are illustrated in the next section by application to the Saint - 
Venant torsion problem. Also shown is how the solution can be extended to elastoplastic 
torsion. The more complicated plane strain and plane stress problems are then con- 
sidered. We conclude with a brief discussion of three-dimensional problems. 

The subject is approached herein from the applied viewpoint and not from the math- 
ematical one. Such questions as the existence and uniqueness of solutions, the differ- 
ences between internal and external problems, problems of convergence, etc. , are not 
discussed. Excellent discussions of many of these questions can be found in references 1 
to 6. 


TORSION PROBLEM 

The torsion problem can be formulated in many ways (refs. 7 and 8): 
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(1) Warping function, W, where 


vV = 0 in R 

( 1 ) 

— = (Zy - mx)a on C 
9n 


Where a is the angle of twist per unit length and Z and m are the direction cosines 
of the outward normal to 'the boundary C of the region R. 

(2) Conjugate of warping function, V, where 


V 2 V = 0 in R 

V = — (x 2 + y 2 ) on C 
2 


( 2 ) 


with G equal to the shear modulus. 

(3) Stress function, F, where 


V 2 F = -2Gce in R 
F = 0 on C 


(3) 


All three formulations represent classical problems in potential theory. This is why 
the methods to be discussed are frequently referred to as potential methods. Number (1) 
is the classical Neumann problem for Laplace’s equation. Number (2) is the Dirichlet 
problem for Laplace’s equation. And number (3) is the Dirichlet problem for Poisson’s 
equation. In what follows, these three equations are referred to as formulations (1), (2), 
and (3). 


Formulation (1) 

As is well known, the solution to Laplace’s equation can be given in terms of either 
a single or double layer potential (refs. 1 and 2). For the warping function we take a 
single -layer potential, that is, 

W(P) = J a (q) log r pq dq (4) 
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where (as shown in fig. 2) P represents an interior point, p and q are boundary points , 
and a (q) is an as yet unknown boundary density function. Interior points are always des- 
ignated by capitol letters, and boundary points by lower case letters. Once a (q) is de- 
termined, the warping function can be immediately calculated at any point. Incidentally, 
the point P can also lie on the boundary. 

In order to determine the density distribution a (q) , use is made of the boundary con- 
dition. Differentiating equation (4) and substituting into the second of equations (1), gives 

W’(p) = f <x(q)|log r 
*/ c 




7r<r (p) = (ly 


mx)pa 


where primes have been used to indicate normal derivatives and (as indicated in fig. 2) 
the normal derivative is taken positive outward. The subscript indicates the point at 
which the normal derivative is calculated. The term no (p) appearing in equation (5) 
arises from the fact that, although W(P) is continuous throughout the whole plane, its 
derivative suffers a jump at the boundary (ref. 1). 

Equation (5) is an integral equation (Fredholm equation of the second kind) for the un- 
known boundary density function o (q). This is an example of an indirect potential method 
since the solution is given in terms of a density function that has no physical significance 
for the problem. 

In order to solve equation (5) for the function a(q), the integral is replaced by a 
summation. The boundary c is divided into m intervals with a nodal point taken at the 
center of each interval. The function o is assumed constant over each interval. Equa- 
tion (5) then becomes 

n 

( A ij - Sy^i = B j 3 ~ l,n (6) 

i=l 

where 



Bj = (l y - mx)j a 

and where r.. is the distance from the midpoint of the j* h interval to a point in the i th 
interval and the derivative is evaluated at the midpoint of the j interval. The integra- 

ii. 

tion is taken over the r interval. 
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The set of n equations is solved for c^. The warping function W can then be com- 
puted at any point by equation (4). We see that there are only n equations to solve for 
the n boundary values of a, rather than for a grid of nodal points throughout the region. 
The integrals over each interval can be evaluated in closed form (for straight boundaries), 
or a Simpson’s rule or even trapezoidal rule may be used. 


Formulation (2) 

Consider now the second formulation. Here we have two alternatives. Either a 
single- or double -layer potential may be used. Since the single -layer potential is 
continuous including the boundary, we have 

V(p) = fa (q) log r pq dq = ^ (xj; + yjj) (7) 

«/ n 


This is a singular Fredholm equation of the first kind for the unknown function a (q). 
Again, the integral is converted to a sum and solved for the at the nodal points. 


where 




i=l 


A ir/ logr ij dq 

+5 f) 


( 8 ) 


Alternately, we can use a double -layer potential and write 

V(P) = 

c 

Since the double-layer potential suffers a jump on reaching the boundary (refs. 1 
and 2), we get 


p(q)(log rpj’ dq 


(9) 
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which is a singular Fredholm equation of the second kind for the unknown density function 
fi( q). Once equation (10) is solved for fi(q), V(P) can be calculated everywhere from 
equation (9). As before, we replace equation (10) by 

v'^Bj (id 

i=l 

where 






Although at first sight it would seem that the formulation in terms of the Fredholm 
equation of the first kind is somewhat simpler, this equation has not been used often. 

The reason is that solutions of Fredholm's first equation are not as well understood as 
those of the second kind. In particular, the equation does not necessarily have a solution 
for every right hand side. The double-layer potential is thus usually used for the 
Dirichlet problem . 

We now come to formulation (3). Here there is introduced a new approach which is 
more general than either the single-layer or double-layer potentials. Furthermore, the 
previous solutions were indirect in the sense that they involved solutions for physically 
meaningless density functions, but the new method is semidirect. 


Formulation (3) 

Equation (3) is an inhomogeneous equation for the stress function F. Neither the 
single- or double -layer potentials will satisfy the inhomogeneous equation because a par- 
ticular solution is needed. However, a solution can be obtained in a different form that 
satisfies both homogeneous and inhomogeneous equations by the use of Green's formula 
or, as it is sometimes called Green's third identity. 
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Given two functions u(x,y) and v(x,y) with continuous first and second derivatives 
in the region R, Green's second theorem states 



(uV 2 v - vV 2 u)dx dy = 




( 12 ) 


If we let u = F and v = log r, this equation results in (see the appendix) 

2jtF(P) = JJ f(Q) log r p Q dx dy + J jr^log rp^’ - F* log r pq Jdq (13) 


where for our problem f(Q) = -2Ga. 

This equation gives the values of F at any interior point p in terms of F and F* 
on the boundary. We notice that the contour integral represents both a single-layer and a 
double-layer potential, but, instead of being artificial density functions, they are propor- 
tional to the distributions of the stress function and its derivative on the boundary. We 
might therefore designate this type of solution as semidirect. 

Since the double layer potential jumps as we approach the boundary, equation (13) be- 
comes for P on the boundary 

^ F (p) = J J f (Q) log FpQ dx dy + J ^(log r^’ - F* log r pq jdq (14) 


This solution has many advantages. 

(1) It works for an inhomogeneous equation. This is particularly important for the 
elastoplastic problem as will shortly be seen. 

(2) It takes care of any boundary condition. If F is prescribed (Dirichlet problem), 
equation (14) becomes an integral equation for F* on the boundary. If F* is prescribed 
(Neumann problem), it becomes an integral equation for the unknown F on the boundary. 
If F is prescribed over part of the boundary and F on the other part (mixed boundary 
value or Hilbert problem), then equation (14) is an integral equation for those parts of F 
and F* that are unknown. 

(3) The unknown functions are physically more meaningful. 

The solution of equation (14) can be obtained exactly as previously described by re- 
placing the integrals by summations. 

To calculate the stresses from the stress function, it is not necessary to carry out 
any numerical differentiation of the stress function. We can differentiate under the inte- 
gral sign in equation (13). Then we can calculate the stresses by the same type of numer- 
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ical integration as before, once F and F' are known on the boundary. This is a more 
accurate way of computing the stresses. 

The use of Green’s formula is of course not restricted to formulation (3) for the 
stress function. It can just as well be used for determining W or V. We just replace 
F by W in equations (13) and (14) and set f equal to zero. This solution is then re- 
ferred to as a direct potential method, since we solve directly for the physical quantity 
desired, the warping function. 


Some Numerical Results 

A number of torsion problems have been solved by Jaswon and Ponter (ref. 9), in- 
cluding prismatic bars with cross sections consisting of solid and hollow ellipses, rec- 
tangles, equilateral triangles, and circles with curved notches. The direct potential 
method was used to determine the warping function. A few of the results taken from this 
reference are shown in table I. Results obtained by the present author for the square 
cross section using the stress function rather than the warping function were essentially 
the same. Excellent agreement is obtained with the known analytical results. 


Elastoplastic Torsion 

By use of Green’s boundary formula we can treat the elastoplastic problem the same 
way as the elastic problem. With plastic flow occurring, the equation for the stress 
function can be written as 


v F = -2Ga - 2G 




zx 


3y 


0X 


(15) 


The function f in the double integral of Green’s formula (eq. (14)) now contains the 
plastic flow terms indicated. These terms are computed by the use of the usual plas- 
ticity theory. The solution is an iterative one. 

We start by assuming the plastic strains to be zero. The problem is then solved us- 
ing Green's formula as indicated previously, and the stresses are computed by differen- 
tiation. One then computes the total strains from the relations 
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and 


(16) 


C = J_ or + gP 

ZX 2 Q zx zx 


1 p 

e = -— a + e 
zy 2G zy zy 


and the equivalent total strains from 


2 J 2 2 

e , = — WeiL + €Z„ 
et r y zx zy 

V3 


(17) 


From the stress -strain curve of the material the equivalent, plastic strain is then 
obtained: 


€ P ~ f ( € et^ 


And the plastic strains are then given by 


and 



(18) 


(19) 


First approximations to the plastic strains and consequently for the function f(x,y) ap- 
pearing in equation (15) have thus been determined. The process is then repeated until 
convergence is obtained. A more detailed discussion of the successive approximation 
method (or the ’’method of initial strains, ” as it is sometimes called) can be found in 
reference 10. Note that only the area integral, which represents the right hand side of 
the set of equations to be solved, changes from iteration to iteration. Note also that, al- 
though total plastic strains have been used in the previous equations for ease of writing, 
incremental strains could have been used just as well. 

To illustrate the correctness of this type of formulation we can consider the problem 
of a circular shaft for which the solution can be obtained in closed form if linear strain- 
hardening of the material is assumed. 
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Consider a circular shaft of radius a. The radial coordinate will be designated by 
p, to distinguish it from r, the distance between the fixed point and the variable point ap- 
pearing in the previous formulas. In polar coordinates, because of symmetry, the func- 
tion f appearing in equation (15) becomes 

f ■ - 2G “ + 7 <20) 

For linear strain hardening, it follows that (ref. 10, p. 255) 

€^0 = Ap + ® (21) 


where 


A = 


3a 


2 


3 + 2(1 + v) 


m 

(1 - m)_ 


-l£ (1 + v)e Q 
B = — 

3 + 2(1 + v) — 

(1-m) 


where v is Poisson’s ratio, m is the strain hardening parameter (ratio of the slope of 
the strain hardening line to the elastic modulus), and €q is the yield strain. 

On the boundary F(a) = 0 and, because of axial symmetry, F*(a) = a constant. 
Green’s boundary formula (eq. (14)) then becomes 


0 = 2G 


// (2A -“ + ! 


log r p Q dx dy - F*(a) j log r pq dq 


which, upon solving for F’(a), gives 


( 22 ) 


F’(a) = G[a(2A - a) + 2B] 


Hence, for any interior point 


2ttF( p ) = 2G A - a + -j log r pq dx dy - F’(a ) £ log r pq dq 


(23) 
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or 


F(p) = | [(2A - a)(p 2 - a 2 ) + 4B (p - a)] 


and the shear stress r is given by 


T = - H = 2G 

Q. 

l 

dp 

L\2 / J 


( 24 ) 


which agrees with the solution obtained in an entirely different fashion in reference 10. 
Note that this solution is valid only in the plastic region, that is, for p > p where p , 
the elastic -plastic boundary, is given by (ref. 10), 


P c 


2(1 + v)e Q 


For p < p c the usual elastic solution prevails. 


(25) 


THE PLANE PROBLEM 


We now consider the plane problem of elastoplasticity, that is, plane strain or plane 
stress. The problem can be formulated in several ways: (1) by the use of the Airy 
stress function, which is a semidirect formulation, (2) by the use of singular solutions of 
the Navier equations, which is a direct formulation, or (3) by the use of fictitious loads, 
which is an indirect formulation. The second approach lends itself directly to the three- 
dimensional extension. The last method cannot easily take into account plastic flow. 


Formulation in Terms of Airy Stress Function 
In terms of the Airy stress function, we have to solve the problem (for plane stress), 

V 4 F = f(x,y) in R (26) 

where F and 8F/3n are given on the boundary, and 
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( 27 ) 


2 l* 2 4 d2e v 

f(x, y) = -Ev(afT) - E[_^ + y 


3y 


3x' 


0 2 € P 

xy 

ax ay 


For plane strain this function is slightly different (ref. 10). For an elastic isother- 
mal problem f(x,y)=0. Let 


<p = v 2 f 


(28) 


Then, 


V 2 <p = f(x,y) 


(29) 


We now have a Poisson equation for cp and an inhomogeneous biharmonic equation for F 
Using the same Green's boundary formula as before for (p results in 


tf<?(p) = JJ l0g r pQ d7? + J [^ 0g r pq)q " * log r pq dq ( 30 > 


For the biharmonic equation we can obtain a similar boundary formula using Green's sec- 
ond theorem (see the appendix) : 


4,rF(P) = fL *> + /H V %q)’ - F ' V Vl + (Ppq)> ' w ]*' 1 m) 

where 


2 , 

p = r log r 

We have thus two coupled integral equations for the two unknown functions (p and 
<p on the boundary. These equations can be solved as before. The boundary is divided 
into m intervals, and F, F*, <p, and <p are assumed constant over each interval. 

The integrals are replaced by summations, giving 
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where 


AU 

n(p i +b i] (p ] ) +J i i = X,2„..,m 
3=1 

m 

4 * F i + d ij^ + e i] F i + f ij F j ) + K i 

3=1 


f 3 log r. . 
a ii = / d( l 

X J 7. Pin 


Cjj = f — d <! 

1 J i 3 “ 

f av2 pij 

e y ,J _i> d, 

3 

J i = ff„ *(Q)logr 1Q d{di) 


b ij = - J. l0 « r lj d « 


d ij = ' jf Pij d< l 


f ii =-/ Aj dq 


K i = ffv, ,(Q)p iQ d? *> 


We thus have 2m equations to solve for the 2m unknowns. If part of the boundary 
happens to be a line of symmetry, then F’ and cp are zero there, and the number of 
unknowns is reduced. Once this solution is obtained, the stress function F at any inte- 
rior point is obtained from (see the appendix) 

8 jtF(P) =J[ p pQ t(|,„)d« dt, + f jV(v 2 ppq)' - F'v 2 p pq + (p pq y v - p Pq d>']dq ( 33 ) 
R * c ^ ^ 

As a simple trivial example, consider an elastic square plate under uniform tension 
with unit load (as shown in fig. 3). By taking just four intervals (each side is one inter- 
val), we get 


cr x = 1.0000 
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We now consider a problem that is by no means trivial, namely, the elastic thermal 
stress in a square plate under a parabolic temperature distribution. The solution has 
been obtained in various ways, the most accurate probably being a finite difference solu- 
tion using a 20 by 20 grid and thus involving the solution of 400 simultaneous equations. 

A comparison of that result (ref. 11) with the present method using 60 intervals on the 
boundary is shown in figure 4. Actually, because of double symmetry, there are only 
15 equations to solve. As shown by the figure, good agreement was obtained. 

A similar formulation has been used for plate bending problems in reference 12. 

Let us now consider a problem for which no good solution has been available here- 
tofore, namely, the elastoplastic problem of an edge -notched beam under pure bending 
(fig. 5). 

The function f(x,y) appearing in equation (31) is the right hand side of equation (27) 
with T = 0. An iterative solution is performed as for the torsion problem except that the 
plane stress or plane strain plasticity equations are used (ref. 10). 

Figure 5 shows the spread of the plastic zone with increasing load for the plane 
strain case. A complete discussion of this problem is presented in references 13 and 14. 


Formulation in Terms of Potential Functions 

The preceding problems can also be formulated in terms of single -layer density 
functions as for the torsion problem. This approach has been investigated in some detail 
in references 3 to 5 and 15 to 17. This formulation is based on the fact that a biharmonic 
function can always be represented in terms of two harmonic functions as shown in stand- 
ard elasticity books. Thus, 


F(P) = r|a(P) + /3(P) 


(34) 


where a and /3 are harmonic functions and can therefore be represented in terms of 
single- or double-layer potentials and rp is the distance from the point P to the or- 
igin of coordinates. We can thus write 


«(P) = M(q) log r Pq dq 

|3(P) = f or(q) log r pq dq 


( 35 ) 


and substituting into equation (34) on the boundary gives 
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> ( 36 ) 


F(p) = Tp J* ju(q) log r pq dq + J a (q) log r pq dq 




F’(p) = -7T Jr p ji(p) + a (p)] + 2r p r p J pi(q) log r pq dq 


+ r* 


J M(q)(log rp^ dq + j a(q)(log r pq ^ dq^ 


which results again in two coupled integral equations for the unknown density functions 
jj. and a . These equations can be solved as previously indicated. A number of solu- 
tions using this approach have been obtained by Jaswon and by Rim as previously noted. 
Alternately, we can use Green's boundary formula for the two potentials. Thus 


ira(p) = 


jr/3(p) = 




a’(q) log r pq dq 


0’(q) log r pq dq 


(37) 


With F and F* given on the boundary, we can write 

|3(q) = F(q) - r q a(q) 

/3’(q) = F’(q) - [r q a(q)] 

Eliminating /3 and /3* from the second of equations (37) gives 


(38) 



r pq) q ' dq - / { P ' (q) - [ r q“ (q> ] } l0g r pq dq (39) 


We thus have, as before, a pair of coupled integral equations for the two unknown func- 
tions a and a , which can be solved as previously . 

A variation of this formulation that is advantageous is given in reference 18. There- 
in a and a are eliminated from equation (39) by using the first of equations (35) and 
its normal derivative on the boundary. This gives an integral equation in just one un- 
known density function ju(q), namely, 
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7T 


F(p) 



m( q) log r p q 




F(q) - Tq 


/ 1o « r qqi 





‘‘(‘•l 1 log r Mi 


dqj 



+ ‘-nr q ii(q) 


log r pq dq 


(40) 


Equation (40) contains only one unknown function: jn(q)„ It can be solved numerically 
(as previously indicated) by dividing the boundary into n segments and writing the inte- 
grals in summation form. This solution yields n simultaneous equations for the n un- 
known values of jz. Once \x is known, the function a is calculated from the first of 
equations (35). 

To determine the function /3, its values on the boundary /3(p) are first calculated 
from the first of equations (38). Then expressing /3 in terms of a double -layer potential 
(similar to eqs. (9) and (10)), we have 


0(P) = 


0(P) - 




where a (p) is an unknown density function, which is now determined by solving the sec- 
ond of equations (41). This again requires a solution of n simultaneous equations as be- 
fore. Once a (p) is known, j3(P) is determined from the first of equations (41). TMs 
completes the solution since F can now be computed everywhere from equation (34). 

The advantage of this approach, even though it is somewhat more complicated, is 
that, instead of having to solve 2n simultaneous equations in the case of n boundary in- 
tervals, we solve separately two sets of n simultaneous equations each. This is advan- 
tageous .from the viewpoints of computer running time and storage, loss of significant 
figures, and round-off errors. 
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Formulation in Terms of Fictitious Loads 



irei, iy). rnis metnoa maKes use oi t 
of the region in such a way that the boundary conditions are satisfied. These fictitious 
lo 


analogous to the potential of a double layer. 

For these fictitious loads one can make use of the fundamental stress state called 
the ’'radial simple stress state” due to a concentrated force P acting on the boundary 
of the half space. If we consider a distribution of these forces— P(q) -acting on the bound 
ary of the body, then the resultant stress vector acting on a plane through an interior 
point M (as shown in fig. 6(a)) will be given by (ref. 19) 


P(q) — ± cos £ e dq 


where (p is the angle between 


e is the unit vector in the direction of r 


and r, a is the angle between n and r and where 
an of r. It should be emphasized that the loads P(q) 




where 


A(P) P x r x + P y r y 

is a function of the fictitious loads P and where 

q -Vx + Vy 

PqR-—-7— 

r 

is a function only of the geometry. Again, the integrals can be replaced by summations. 

A similar type of formulation was presented by Liu and coworkers (refs. 20 to 22), 
who also considered multiply connected regions, bodies with axial symmetry, and ex- 
tended the formulation to three-dimensional problems as well. We note that the method 
is not directly applicable to bodies with displacement boundary conditions. A similar 
type of approach has also been used by Oliveira (ref. 23) using complex potentials. 


Direct Formulation 


We now come to a direct formulation of the two -dimensional problem which is 
directly extendable to three -dimensional problems. This formulation is based on a sin- 
gular solution of the Navier equations. 

The Navier equations for plane strain, including plastic flow, temperature , and body 
force terms, are 


and 


V 2 u. + — - — 6 2 = 2ep. . + 2 1 + - (ofT) . - 
1 l - 2v ’ 1 - 2v > l G 


0H£ ]3 = u ],3 M = 1 ’ 2 


(45) 


Henceforth the usual tensor notation will be used, where a repeated subscript indicates 
summation over its range and a comma indicates partial differentiation. In equation (45) 
€?• represents the sum of all the plastic strain increments up to and including the cur- 
rent increment of load, v is Poisson’s ratio, a is the coefficient of linear thermal ex- 
pansion, and G is the shear modulus. 
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For plane stress the equation becomes 


vV + I±± e 


= 2ef, , + ^>L e P. + W ..! 1 

13 > 3 1 _ 7. >* 1 _ V > X G 


(45a) 


where 


eP = + £ p 

17 x y 


The boundary conditions that must be satisfied by the displacements over that part of 
the boundary where forces P. are specified are 


G 


= (u, 4 + u. -)n. + — u k k n. - 2 (e?.n. + -—t i l aTnA 

hl 3 1 -2v K,K \ 13 3 1 - 2v / 


for plane strain and 


p i 


= (a i,J + “m'”) + 777"k,k n l 


to 


-0—e p . i±± «T)n, 
A - v 1 - v ) 1 


(46a) 


for plane stress. 

We now note the following interesting fact: By defining a pseudobody force by 


F. = F. - 2G 


e l1 1+ i 

13 ’ 3 1 - 2v > l 


(47) 


and a pseudoboundary force by 


P.’ = P. + 2G(e?.n. + 1+ v aTn.] 

1 1 V 133 1-2* l ) 


(47a) 


we may write equations (45) and (46) for plane strain as 


2 i f: 

vV + — L_ e=.i 

1 1 - 2* ’ l G 


> 


P 


2* 

G '~ 1 ’ 3 ’ “ 3 ’ r "’ 3 ' 1 - 2* 


- = K , + u. 4 )n. + 


Vk ! 


n. 


(48) 
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We thus have to solve a problem with modified body forces and modified boundary 
loads. This is a direct generalization of the Duham el -Neumann analogy for thermoelas- 
tic problems (ref. 24),, Equations (48) could have been obtained by assuming an elastic 
body with body forces given by equation (47) and surface forces given by (47a) (for plane 
strain). For this imaginary body the equilibrium equations would contain the plastic flow 
terms together with the body forces, if any, but the stress -strain relations would not. 
The displacements obtained by solving equations (48) would thus be correct. The 
stresses, however, would be pseudostresses. The true stresses have to be computed 
from the correct stress -strain relations. It should also be emphasized that the actual 
pseudobody and pseudoboundary forces are not known apriori and have to be determined 
as part of the solution by, for example, an iterative process. For plane stress, P f and 
F* become 


F = F. - 2G 


P = P. + 2G 

I 


cf ; + _iL. e p + l±±(oT) i 

X J5 J 1 _ V > 1 1 _ V ’ X 

Sn, + (-J^s p -!±i:tfrk 


A 




'll 3 


,1-y l-i/ 


(49) 


A solution of the nonhomogtneous Navier equations can be obtained by making use of 
the singular solution of the Navier equations due to a point load, as given for example in 
Love’s classical book on elasticity (ref. 25), and also making use of Betti’ s reciprocal 
theorem. We then arrive at a solution (see the appendix) that is known as Somigliana’s 
identity, namely, 


XU 1 (P) - / (UyPj - TyUjWq ♦ 


U..F. dA 
R *3 3 


(50) 


where 


6ji(3 - v) log r - r .r , 

_ y 9 * ? 3 j 


T.. = - 


13 477 (1 - v)r 


U.. = - x 
13 8 ttG(1 - v) 

1 {[yi - 2v) + 21 -jrJ |-C1- 2,Kr ;1 nj - r ,»,)} 




( 51 ) 


and F ? and P* are given by equations (47) and (47a). The coefficient X is equal to 1, 
if P is an interior point and is equal to 1/2, if P = p is a boundary point. For plane 
stress, v in equation. (51) is replaced by v/ ( 1 + v), and P ? and F* are taken from 
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equations (49), The derivatives appearing in equation (51) are taken at the variable point 
of integration. 

Instead of using the pseudobody and pseudoboundary forces in equation (50), an alter- 
nate solution can be written as follows (see the appendix) 

Au.(P) = £ (U ijPj - T..u.)dq + U..F. dA + S jki (eP k + S^o^dA (52) 


where 


s jki ~ 


4tt(1 - v)r 


[( 6 ij r ,k + 6 ik r ,j " 6 jk r ,i )(1 ' 2v) + 2r ,i r ,j r ,k] 


(53) 


It can be shown that equations (50) and (52) are identical. Formulation (52) has the 
advantage that it uses the physically given boundary conditions and, more importantly, 
that it does not require the derivatives of the plastic strains. 

The solution is now obtained by replacing the integrals by sums as before. Note that 
the displacements and, hence, the stress fields are given directly in terms of the bound- 
ary forces and boundary displacements. If the displacements are given on the boundary, 
we have an integral equation for the boundary forces. If the forces are given on the 
boundary , we have an integral equation for the unknown boundary displacements. If the 
forces are given over part of the boundary and the displacements over the rest of the 
boundary (mixed boundary value problem), we have an integral equation for those parts of 
the boundary loads and displacements that are unknown. Also, since we are working di- 
rectly with displacements, there is no difficulty with multiply connected regions. 

The stresses at any interior point P are obtained by direct differentiation of equa- 
tion (52), resulting in 

dA 

(54) 

where 






' ijk u k 


R 


V... T F, 

i]k k 


dA 


n 


R 


T f 
ijkZ 1 


~kZ 


r) 


21 



( 55 ) 


-\ 


V ijk ” G (yik,j + U jk,i + 5 ij t _ 2 j/ U Zk,Z 


T i ik =- G L^ M + 5 ^\ > 


T» — — G 1 Si, i • • + Zi_i * • + 5|* ' " Si-y™ rvi 

ijkZ \ ktijj kZ j , i lj 1 — 2js kZiUjin 


As previously noted all derivatives are taken at the variable point of integration. For the 
elastic problem all the plastic flow terms disappear. 

Solution to a number of elastic problems using these equations have been obtained, 
for example, in references 6 and 26 to 30. These include problems of circular and ellip- 
tic regions, rings, rectangles, bodies with holes, inclusions, notches, and cracks as 
well as some thermoelastic problems. 

The Navier equations (45) and (46) and the corresponding solution (52) have been 
written in terms of displacements, stresses, and plastic strains. This is the appro- 
priate form if the method of successive elastic solutions, or method of initial strains 
(ref. 10), is used to solve the elastoplastic problem. If the tangent modulus method is 
used (refs. 31 to 33), these equations are written in terms of velocities, stress rates, 
and plastic strain rates. The boundary -integral formulation remains the same. 


THREE-DIMENSIONAL PROBLEMS 

For three-dimensional problems the extension is direct. The Navier equations (45) 
retain the same form, except that the range of subscripts is now three and the Laplacian 
appearing in these equations is the three-dimensional Laplacian. The singular solution, 
instead of having a logarithmic singularity, has a 1/r singularity. The solution then 
takes the same form as before, except that the kernels of the integrals are slightly dif- 
ferent. The area integral becomes a volume integral, and the line integrals become sur- 
face integrals over the bounding surface of the body. Thus equation (52) becomes 

XUi(P) - f s (UyPj - VjldS + f v UjjFj dV Sjkl ( e j> k ♦ V^dV (56) 


where 
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u., + - i [(3 - 4 v) 6-j + r .r .1 

13 1677(1 - t*)G r L 13 ’ l ’ 3J 


A 


T-- = - - 

13 877(1 - v) r 2 


- L/[(l . 2 v) 6.. + 3r r .1 — + (1 - 2v)(r n. - r n,)| l 

. - v ) 2 l L 13 - 1 ’ 3 -l 0n > ] 1 > x 3 J I 


2 .. . = 

jki 


1— ± [(1 - 2,)( 6ij r >k + V,j - V,i' + 3r t l r ,i r ,k] 


877(1 - v) r 2 


( 57 ) 


Instead of dividing the boundary into a set of linear intervals, we now have to divide 
the bounding surface into a series of surface elements. These may be rectangular or 
triangular in shape as in the two-dimensional finite -element formulation. The unknown 
functions are assumed constant over each element, and the integrals are replaced by 
sums over all the elements. If the surface is divided into m elements, we get 3m 
simultaneous equations to solve. 

The stresses are given by 


<r ii (p) = I (v iik p k - + L v ijk F k dv - 2& ij - L T m (4i + 


(58) 


where T.^, T.^^ are as given by equations (55) and (57). The general elasto- 

plastic flow theory leading to these equations is given in reference 34. 

Using this approach solutions have been obtained for several elastic problems in ref- 
erences 27 and 35 to 37 among others. Figures 7 and 8 show two such problems. Fig- 
ure 7 (from ref. 35) shows the results for the simple problem of a cube loaded in tension 
at one end, while on the other end, as well as on two normal faces, the normal displace- 
ment components are set equal to zero. The results are shown for 12 surface elements, 
each surface being divided into two triangles. Good agreement is obtained with the exact 
solution. 

An example of a much more difficult problem involving a semielliptical surface flaw 
in a plate subjected to tension is shown in figure 8 (from ref. 37). An exact solution is 
not available here for comparison. A similar problem is treated in reference 30. 


CONCLUDING REMARKS 

This brief survey indicates the potential of the boundary integral methods for solving 
two- and three-dimensional elastic and elastoplastic problems. The more important ad- 
vantages of these methods are the reduction of the problem size, since the unknown quan- 
tities appear only on the boundaries, and the ability to handle arbitrarily shaped 
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boundaries. This advantage is particularly apparent in using the direct -method involving 
the Navier equations, since neither smooth boundaries nor simply connected regions are 
required. 

The methods described are not limited to elastostatic or elastoplastic analysis of 
isotropic materials. References 38 to 40 formulate the general transient elastodynamic 
problem. Reference 41 discusses the linear viscoelastic problem; references 42 and 43 
present solutions to anisotropic elastic problems. Problems with inhomogeneities are 
discussed in references 27 and 44. 

Finally it should be pointed out that, although the theory pertaining to these methods 
has been explored in some detail, the numerical applications are still in their early 
stages. Thus, for example, reference 29 discusses improvements that might be made by 
assuming linear variations of the unknown functions over the boundary intervals rather 
than assuming them constant over each interval. Reference 13 indicates that taking some 
intervals to be gradually varying in size, rather than having them all the same size, can 
improve both the stability and accuracy of the solution. The question of using iterative 
solutions of the integral equations involved (as was done, for example in ref. 19), instead 
of reducing them to sets of algebraic equations, has not been fully explored. Thus a 
great deal of work remains to be done in developing the proper numerical techniques for 
using these powerful boundary integral methods to their utmost potential. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, July 24, 1973, 

501-21. 
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APPENDIX - SOME MATHEMATICAL DERIVATIONS 


Derivation of Equation (13) 


Substituting u = F and v = log r into equation (12) gives, since v = 0 and 
V 2 F = -2Ga = f(Q), 



f(Q) log r p Q dx dy = 



F’ log r pq dq 


(Al) 


Now the integrands in equation (Al) are singular when r = 0. We must therefore ex- 
clude this point from the region R by drawing a small circle of radius e about this 
point and making a cut as shown in figure 9. The line integration is performed as shown, 
and the limit taken as e shrinks to zero. Thus, 


and 




27rF(P) as e —0 


(A2) 



/ 


log r — dq - e log e 
9n 



as e — 0 

(A3) 


which results in equation (13). 

For a point on the boundary a semicircle must be used instead of a circle to exclude 
the singular point leading to a 7rF(p) term instead of a 27 tF(P) term and resulting in 
equation (14). 
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Derivation of Equations (31) and (33) 


2 

Substituting u = V F ir«.o equation (12) results in 


// V 2 F V 2 v dx dy = / 
R 


vV F dx dy + 


/[ v 


2 F - V 9(V F) 
9n 8n 


dq (A4) 


Since the left side is 'symmetric in F and v, we can also write 




V 2 FV 2 v dx dy 


2n 


= [[ FV 4 v dx dy + /|>vH-fMTx) 
yy R L 3n 3n 


dq (A5) 


Subtracting the second equation from the first gives 


IL 


4 4 

(FV V - VV 


F)dx dy = /‘[Fi<5!v).^V 2 v-vlff!e + Slv 2 Fldq (A6) 
7 L 3n 3n 3n 0n J 


9 4 

Let v = r log = p where r is as previously defined. Then V v = 0, and, since 

V 4 F = f(x,y), 



f(£,77)p(x,y,£,r?)d£ d ?7 = 



9(V 2 p) _ 3F 
0n 9n 



P 


g(V F) + 0p v 2 f 


3n 


3n 


dq (A7) 


As before, the singular point in the region must be excluded before evaluating the line 
integral. This is done by drawing a small circle of radius e about this point, making a 
cut (as shown in fig. 9), performing the integration shown, and taking the limit as e 
shrinks to zero. Only the first term of the line integral makes a contribution. Thus 

F i e de = -8 ttF(P) (A8) 

€ 

Equation (A7) then becomes 


lim 

€-0 


/ 


F d( V £) dq = 
9n 


= -lim f 
e-0/ 
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8irF(x,y) = 


Ik 


pf(^,7j)d| drj + f F MI p) - v 2 p + V 2 F 
r /, L 0n 0n 8n 


0(V 2 F) 
0n . 


dq (A9) 


For a point on the boundary a semicircle is drawn to exclude the singular point re- 
sulting in the term 47rF(p) rather than 8 ttF(P). 


Derivation of Equations (50) and (5 1) 


Betti’s reciprocal theorem states: Given two states of stress and displacement due 
to two different sets of loads, then 


/ P,u* dq + / F.u* dA = / P,*u, dq + / F.*u. dA 
J c J J Jji i J J c J j J J 


(A 10) 


For the problem under consideration, F. and P. are taken to be equal to F.’ and P.' 

j j 3 3 

(as given by eqs. (47) or (49)), and Uj is the displacement field for the body. For the 

second stress state we choose Kelvin’s singular solution for a point load in an infinite 
medium. This solution (refs. 6 and 25) is given by equations (51) where U^. is the dis- 
placement in the Xj direction due to a unit load in the x^ direction a distance r away, 
and Ty is the corresponding traction tensor at the same point. Thus, for a unit con- 
centrated force acting in the direction given by the unit base vector e^, 


and 


* 

u . = U, -e, 

] kj k 


P j =T kj e k 


(All) 



= 0 J 


Substituting into equation (A10) gives 

4 v d? + 4 e ^ = 4 v dq <Ai2 

where R £ with boundary c £ is the circular region excluded because of the singular na- 
ture of the Kelvin solution. Taking the limit as e goes to zero, in a manner similar to 
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what was done in deriving equation ( 13) , results in equation (50). The term Xu.(P) 
comes from the term on the right hand side of equation (A 12). 

Equation (52) is obtained as follows. The elastic part of the total strains is given by 


Now 


e® = e.. - eS - 6.. aT 
U i] i] i] 


/ 


* e 

dA 

R.R e 13 13 


■I. 


dA 

R-R, 13 13 


(A 13) 


(A 14) 


where the starred fields refer, as before, to Kelvin's singular solution. Equation (A14) 
follows from Hooke’s law. Substituting equation (A13) into equation (A 14) and using the 
strain displacement relations, the equilibrium equations for the starred stress field, and 
the divergence theorem result in 


/ ¥ii”i d( l - / ff iK € S + 6 ii = / u i d q + / u.* F. dA 

*x+c _ 1 13 3 ^R-R, 13 v 13 13 ' *x+c 1 13 3 •'R-R 1 1 . 


Utilizing the relations 


a i*j n j = P i* = T ki e k 


u i = JJ ki e k > 


a.. = S. M e 1 
ij ijk k 


(A 15) 


(A 16) 


where T^, u^, and are defined in equations (51) and (53), gives 

/ w*- / R v( e ii +5 ii^h=£ v,-dA 


(A17) 


Taking the limit as e approaches zero results in equation (52). 
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Figure 3. - Square plate under uniform tension. 
Exact solution: o x = (p - 1. Boundary-integral 
solution, o x = 1. 000. 



Figure 4. - Thermal stress in square plate; Oy at y = 0 
as function of x. (One quadrant of plate shown. ) 
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i 5. - Growth of the plastic zone size with load for a linear strain hardening specimen with a 10° 
notch subjected to pure bending; plane strain; ratio of slope of strain hardening line to elastic 
, 0.05 (from ref. 13). 





P(q) 

(a) Stress vector at interior point M 
due to fictitious load distribution. 


y 



(b) Boundary forces and variables entering into 
fictitious load formulation. 

Figure 6. - Geometric and force variables enter- 
ing into fictitious load formulation (eqs. (43) 
and (44)). 
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Figure 7. - Unit cube under axial load 
for 12 surface elements. 
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(a) Geometry and relative dimensions. 

Figure 8. - Semielliptical surface flow in a plate subjected to tension (from ref. 37). 


Distance from plate center, xj 

(b) Axial stress distribution along major axis. 256 Surface elements 
used. 



Figure 9. - Integration path excluding singular point. 
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